# build_parameters.R
#
#

# ------------------------------------------------------------------------------
# set general parameters
# ------------------------------------------------------------------------------

# relative path to where the data files are stored
data.dir = file.path("build","data")

# name of file to write the elasticities
output.file = file.path("data","parameters.gms")

# when aeo data is used these define the aeo release year and scenario to use
aeo.year = 2018
aeo.scenario = "REF2018"

# years to include in the model
sage.years = seq(2016,2061,by=5)



# ------------------------------------------------------------------------------
# load libraries
# ------------------------------------------------------------------------------

library(httr)
library(RJSONIO)

# load in the R utilities for sage
source(file.path("utilities","R_utilities.R"))



# ------------------------------------------------------------------------------
# armington elasticities
# ------------------------------------------------------------------------------

# load the gtap elasticities
gtap = read.csv(file.path(data.dir,"hertel et al - GTAP v7 - behavioral parameters - table 14-2.csv"),
                stringsAsFactors=FALSE,header=TRUE,skip=21)

# load the gtap to sage mapping
gtap.map = read.csv(file.path(data.dir,"gtap_to_sage.csv"),
                    stringsAsFactors=FALSE,header=TRUE,skip=1)

# get the list of sage sectors
sage.sectors = unique(gtap.map$sage)

# for each sage sector get the armington elasticities
se_nf = NULL
se_dn = NULL
for (s in sage.sectors) {

  # gtap sectors within the current sage sector
  ind = gtap$GTAP.Sector %in% gtap.map$gtap[gtap.map$sage==s]

  # compute the value weighted national-foreign armington elasticity
  # weighting is by import value at world prices as reported in GTAP v9
  se_nf_value = sum(gtap$sigma_D[ind]*gtap$viws[ind])/sum(gtap$viws[ind])

  # crude oil and natural gas elasticities from hertel et al. (2007) are less
  # appropriate for the U.S. and therefore, we follow caron and rausch (2013)
  # in there updating of the hertel et al. (2007) elasticities for these sectors
  if (s=="cru")
    se_nf_value = 7.3
  if (s=="gas")
    se_nf_value = 2.8

  se_nf = rbind(se_nf,cbind(s,se_nf_value))

  # coughlin and novy (2013) ratio of intra-national to international border effects
  delta = 1.868

  # caron and rausch (2013) relationship between domestic-national and national-
  # foreign armington elasticities
  se_dn = rbind(se_dn,cbind(s,.5-delta*(.5-se_nf_value)))

  # based on Ross (2014) DIEM domestic-national armington elasticities
  #se_dn = rbind(se_dn,cbind(s,3))

}



# ------------------------------------------------------------------------------
# transformation elasticity of output
# ------------------------------------------------------------------------------

# domestic vs. export transformation elasticity of output
# taken from us-rep
te_dx = 2;



# ------------------------------------------------------------------------------
# transformation elasticity for sector differentiated extant capital
# ------------------------------------------------------------------------------

# transformation elasticity for sector differentiated extant capital
te_k_ex = 1.5;



# ------------------------------------------------------------------------------
# klem substitution elasticities
# ------------------------------------------------------------------------------

# load the koesler and schymura (2015) elasticities
klem = read.csv(file.path(data.dir,
                          "koesler and schymura - 2015 - klem elasticities.csv"),
                stringsAsFactors=FALSE,header=TRUE,skip=8)

# for sectors without kl elasticities (either due to inf answers or convergencece
# failure na) and with values in young (2013) use the gmm non-normalized values
# from that study instead
klem$se_kl[klem$nace=="23"] = 0.732      # petroleum refining
klem$se_kl[klem$nace=="19"] = 1.414      # leather goods
klem$se_kl[klem$nace=="E"]  = 1.002      # water and electric utilities

# for sectors with an inf kl elasticity that don't have a direct young (2013)
# correspondence use the services sector value from young (2013)
# value
klem$se_kl[klem$nace=="M"]  = 0.999      # education
klem$se_kl[klem$nace=="50"] = 0.999      # sale and repair of vehicles and fuel

# petroleum refining has an inf value for the kle elasticity, instead use the
# the total industries value
klem$se_kle[klem$nace=="23"] = klem$se_kle[klem$nace=="TOT"]

# load the nace code to sage mapping
nace.map = read.csv(file.path(data.dir,"nace_to_sage.csv"),
                    stringsAsFactors=FALSE,header=TRUE)

# get the list of sage sectors
sage.sectors = unique(gtap.map$sage)

# for each sage sector get the value added elasticity
se_kl = NULL
se_kle = NULL
se_klem = NULL
for (s in sage.sectors) {

  # nace sectors within the current sage sector
  ind = klem[,1] %in% nace.map$nace[nace.map$sage==s]

  # compute the value weighted value added elasticity
  # weighting is by output in the last year of the koesler and schymura sample
  y = sum(as.double(klem$output[ind]))
  kl_value = sum(as.double(klem$se_kl[ind])*as.double(klem$output[ind]))/y
  kle_value = sum(as.double(klem$se_kle[ind])*as.double(klem$output[ind]))/y
  klem_value = sum(as.double(klem$se_klem[ind])*as.double(klem$output[ind]))/y

  # write the value added elasticity to the output file
  se_kl = rbind(se_kl,cbind(s,kl_value))
  se_kle = rbind(se_kle,cbind(s,kle_value))
  se_klem = rbind(se_klem,cbind(s,klem_value))

}



# ------------------------------------------------------------------------------
# interfuel energy substitution elasticities
# ------------------------------------------------------------------------------

# weighted average of primary morishima elasticities from serletis et al (2010a)
# where the weights are the input value for the last year in their dataset
# transportation  based on serletis et al (2010b) international study
#             industrial commercial household electricity transportation
se_en_agg = c(0.33,      0.24,      0.24,     0.30,       0.25)

# weighted average of primary morishima elasticities from serletis et al (2010)
# where the weights are the input value for the last year in their dataset
# transportation  based on serletis et al (2010b) international study
#              industrial commercial household electricity transportation
se_ene_agg = c(0.61,      0.53,      0.67,     0.01,       0.25)

# commercial elasticities apply to the service sectors
com = c("hlt","srv")

# transportation sectors
trn = c("trn","ttn")

# for each sage sector get the energy elasticities
se_en = NULL
se_ene = NULL
for (s in sage.sectors) {

  # use the commercial elasticity index for designated commercial sectors
  if (s %in% com)
    ind = 2
  else if (s == "ele")
    ind = 4
  else if (s %in% trn)
    ind = 5
  else
    ind = 1

  # save sector specific elasticities
  se_en = rbind(se_en,cbind(s,se_en_agg[ind]))
  se_ene = rbind(se_ene,cbind(s,se_ene_agg[ind]))

}

# save the energy elasticities for households
se_cen = se_en_agg[3]
se_cene = se_ene_agg[3]



# ------------------------------------------------------------------------------
# labor supply elasticities
# ------------------------------------------------------------------------------

# labor supple elasticity - substitution effect
# based on most recent cbo literature review midpoints
lse_sub = 0.20

# labor supple elasticity - income effect
# based on most recent cbo literature review midpoints
lse_inc = -0.05



# ------------------------------------------------------------------------------
# household consumption elasticities
# ------------------------------------------------------------------------------

# substitution elasticity across consumption goods and transportation
# from usrep
se_c = 1.001;

# substitution elasticity across non-energy and energy goods
# from usrep
se_cme = 0.25;

# substituion elasticity across non-energy and non-transportation goods
se_cm = 0.25;

# inverse elasticity of substitution of consumption over time
eta = 1.001;



# ------------------------------------------------------------------------------
# fixed factor shares in resource and agriculture sectors
# ------------------------------------------------------------------------------

# fixed factor resource share of capital for resource extraction sectors
# based on EMPAX
res_share = 0.25

# fixed factor resource share of capital for agriculture and forestry sectors
# based on EMPAX
land_share = 0.40



# ------------------------------------------------------------------------------
# resource extraction and agricutlure sector supply elasticities
# ------------------------------------------------------------------------------

# natural gas parameters based on Arora and Cai (2014)

# short-run supply elasticity for resoursce extraction sectors
se_res_sr = cbind(c("col","cru","gas","min","agf"),c(0.40,0.05,0.02,5,1.19))

# long-run supply elasticity for resoursce extraction sectors
se_res_lr = cbind(c("col","cru","gas","min","agf"),c(1.90,.25,0.50,5,1.19))

# rate of convergence to long-run supply elasticity
rate_sr_lr = cbind(c("col","cru","gas","min","agf"),c(0.07,0.15,0.15,0.07,0.07))



# ------------------------------------------------------------------------------
# macroeconomic variables
# ------------------------------------------------------------------------------

# benchmark population
# Macroeconomic Indicators : Key Labor Indicators : Labor Force, Reference, AEO
n = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                   ".ECI_LBF_NA_NA_NA_NA_NA_MILL.A"))

# rate - population growth rate
# based on average labor force growth rate in AEO
gamma = mean(n$value[-nrow(n)]/n$value[-1]-1)

# rate - labor productivity growth
# Macroeconomic Indicators : Key Labor Indicators : Nonfarm Labor Productivity, Reference, AEO
temp = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                      ".ECI_LP_NOFM_NA_NA_NA_NA_Y09EQ1.A"))
omega = mean(temp$value[-nrow(temp)]/temp$value[-1]-1)

# rate - depreciation rate
delta = 0.05

# rate - benchmark interest rate
rbar = 0.045

# if this is a static model run then redefine the rate parameters - in the
# static run these parameters essentially lose their meaning so these values are
# just to ensure the solution represents a static run and the values should not
# be interpreted with the parameters normal meaning for the dynamic model
if (length(sage.years)==1) {
  gamma = 0.0
  delta = 1.0
  omega = 0.0
  rbar = 0.0
}



# ------------------------------------------------------------------------------
# average energy intensity growth rate
# ------------------------------------------------------------------------------

# baseline improvements in the energy intensity of production are based on the
# unit energy consumption in eia's aeo forecast. for each sage sector we compute
# the gorwth rate of consumption per unit of real output aggregated across the
# nems sectors within that sage sector over the aeo time horizon. the mapping of
# nems to sage sectors is located in build/data/nems_to_sage.csv

# sectors represented in nems
nems.sectors = c("ref","fdp","ppm","bch","ggr","cem","ism","aap","fbp",
                 "mchi","cmpr","teq","eei","wdp","pli","bmf","agg","cns","ming",
                 "comm")

# initialize space for the value of shipements and energy consumption data
output = NULL
energy = NULL

# for each nems sector retrieve the data using the eia api
for (i in nems.sectors) {

  # eia is not consistent in its structuring of api series ids across sectors so
  # we need to handle cases that are different from the norm
  if (i %in% c("wdp","pli","bmf")) {
    id = paste0("ECI_SHP_IDAL_",toupper(i),"_NA_NA_NA_BLNY09DLR")
  } else if (i=="cmpr") {
    id = "ECI_VOS_IDAL_CMPEL_NA_NA_NA_BLNY09DLR"
  } else if (i=="fbp") {
    id = "ECI_VOS_IDAL_FABM_NA_NA_NA_BLNY09DLR"
  } else if (i=="comm") {
    id = "ECI_VOS_NIND_NA_NA_NA_NA_BLNY09DLR"
  } else {
    id = paste0("ECI_VOS_IDAL_",toupper(i),"_NA_NA_NA_BLNY09DLR") }

  # add the output data for the sector to the dataset
  temp = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".",id,".A"))
  temp$sector = i
  output = rbind(output,temp)

  # eia is not consistent in its structuring of api series ids across sectors so
  # we need to handle cases that are different from the norm
  if (i %in% c("ism","aap","ref")) {
    id = paste0("CNSM_NA_IDAL_",toupper(i),"_TOT_NA_NA_TRLBTU")
  } else if (i=="comm") {
    id = "CNSM_ENU_COMM_NA_NG_NA_NA_QBTU"
  } else {
    id = paste0("CNSM_NA_IDAL_",toupper(i),"_NA_NA_NA_TRLBTU") }

  # add the energy consumption data for the sector to the dataset
  temp = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".",id,".A"))
  temp$sector = i
  energy = rbind(energy,temp)

}

# load the nems to sage mapping
map = read.csv(file.path(data.dir,"nems_to_sage.csv"),
               stringsAsFactors=FALSE,header=TRUE)

# get the list of sage sectors
sage.sectors = unique(map$sage)

# for each sage sector get the energy intensities
ene_growth = NULL
for (s in sage.sectors) {

  if (any(map$nems[map$sage==s]==""))
    next

  # get energy and output for the nems sectors in the sage sector
  sector.energy = energy[energy$sector %in% map$nems[map$sage==s],]
  sector.output = output[output$sector %in% map$nems[map$sage==s],]

  # aggregate across nems sectors by year
  sector.energy = aggregate(sector.energy$value,by=list(sector.energy$year),sum)
  sector.output = aggregate(sector.output$value,by=list(sector.output$year),sum)

  names(sector.energy) = c("year","value")
  names(sector.output) = c("year","value")

  # average energy intensity growth rate
  intensity = sector.energy
  intensity$value = sector.energy$value/sector.output$value
  growth = log(intensity$value[nrow(intensity)]/intensity$value[1])/(nrow(intensity)-1)
  ene_growth = rbind(ene_growth,data.frame(ss=c("gas","ele","col","ref"),s=s,value=growth))

}

# energy intensity growth for the transportation sector is based on fuel
# efficiency improvements forecast for airline travel as reported in aeo
temp = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".CNSM_NA_TRN_AIR_USE_NA_NA_QBTU.A"))
temp$value = temp$value/get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                                       ".KEI_TRV_TRN_NA_AIR_NA_NA_BLNSEATMLS.A"))$value
ene_growth = rbind(ene_growth,
                   data.frame(ss=c("gas","ele","col","ref"),s="trn",
                              value=log(temp$value[1]/temp$value[nrow(temp)])/(nrow(temp)-1)))

# energy intensity growth for the truck transportation sector is based on fuel
# efficiency improvements forecast for frieght truck travel as reported in aeo
temp = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".CNSM_NA_TRN_FGHT_USE_NA_NA_QBTU.A"))
temp$value = temp$value/get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                                       ".KEI_TRV_TRN_NA_FGHT_NA_NA_BLNVEHMLS.A"))$value
ene_growth = rbind(ene_growth,
                   data.frame(ss=c("gas","ele","col","ref"),s="ttn",
                              value=log(temp$value[1]/temp$value[nrow(temp)])/(nrow(temp)-1)))



# ------------------------------------------------------------------------------
# average growth rate of coal fired electricity generation share
# ------------------------------------------------------------------------------

# electricity output valued at initial year prices
ele = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                     ".GEN_NA_ELEP_CHP_TTG_NA_USA_BLNKWH.A"))
ele$value = ele$value*get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                                     ".PRCE_NA_ELEP_NA_EDU_NA_USA_Y13CNTPKWH.A")
                              ,year=min(sage.years))$value
ele = ele[sort(ele$year,index.return=TRUE)$ix,]

# coal input valued at initial year prices
col = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                     ".CNSM_ENU_ELEP_NA_STC_NA_NA_QBTU.A"))
col$value = col$value*get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                                     ".PRCE_REAL_ELEP_NA_STC_NA_NA_Y13DLRPMMBTU.A")
                              ,year=min(sage.years))$value
col = col[sort(col$year,index.return=TRUE)$ix,]

# coal fired electricity generation share
share = col$value/ele$value

# average growth rate of coal fired electricity generation share
col_ele_growth = log(share[length(share)]/share[1])/(length(share)-1)



# ------------------------------------------------------------------------------
# household energy intensity growth
# ------------------------------------------------------------------------------

# consumption forecast in aeo
c = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                   ".ECI_NA_NA_NA_GDP_RCNSP_NA_BLNY09DLR.A"))

# consumption electricity intensity growth is based on residential electricity
# consumption growth in aeo relative to consumption growth
e = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".CNSM_NA_RESD_NA_ELS_NA_USA_BLNKWH.A"))
c.share = e$value/c$value
growth = log(c.share[1]/c.share[length(c.share)])/(length(c.share)-1)
cd_ene_growth = data.frame(s="ele",value=growth)

# consumption natural gas intensity growth is based on residential natural gas
# consumption growth in aeo relative to consumption growth
e = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".CNSM_ENU_RESD_NA_NG_NA_NA_QBTU.A"))
c.share = e$value/c$value
growth = log(c.share[1]/c.share[length(c.share)])/(length(c.share)-1)
cd_ene_growth = rbind(cd_ene_growth,data.frame(s="gas",value=growth))

# consumption refined oil intensity growth is based on light duty vehicle energy
# intensity growth in aeo
e = get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,".CNSM_NA_TRN_LDV_USE_NA_NA_QBTU.A"))
c.share = e$value/get.eia(paste0("AEO.",aeo.year,".",aeo.scenario,
                                 ".KEI_TRV_TRN_NA_LDV_NA_NA_BLNVEHMLS.A"))$value
growth = log(c.share[1]/c.share[length(c.share)])/(length(c.share)-1)
cd_ene_growth = rbind(cd_ene_growth,data.frame(s="ref",value=growth))

# zero out any change in coal
cd_ene_growth = rbind(cd_ene_growth,data.frame(s="col",value=0))



# ------------------------------------------------------------------------------
# write parameters to the input file
# ------------------------------------------------------------------------------

# open a connection to the output file
fc = file(output.file,"w")

# write the set defining the model's years
writeLines("sets",con=fc)
t = sage.years[1]
if (length(sage.years)>1)
  for (i in 2:length(sage.years))
    t = paste0(t,", ",sage.years[i])
writeLines(paste0("  t                     time periods            / ",t," /;"),con=fc)
writeLines("",con=fc)

# armington elasticities
write.gams.variable(fc,"se_nf",se_nf,comments="substitution elasticity - national-foreign armington elasticity",write.zeros=TRUE)
write.gams.variable(fc,"se_dn",se_dn,comments="substitution elasticity - domestic-national armington elasticity",write.zeros=TRUE)

# transformation elasticity - domestic and exported production
write.gams.variable(fc,"te_dx(s)",te_dx,comments="transformation elasticity - domestic and exported production",write.zeros=TRUE)

# transformation elasticity for sector differentiated extant capital
write.gams.variable(fc,"te_k_ex",te_k_ex,comments="transformation elasticity for sector differentiated extant capital",write.zeros=TRUE)

# klem substitution elasticities
write.gams.variable(fc,"se_kl",se_kl,comments="substitution elasticity - value added",write.zeros=TRUE)
write.gams.variable(fc,"se_kle",se_kle,comments="substitution elasticity - value-added-energy bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_klem",se_klem,comments="substitution elasticity - materials and value-added-energy bundle",write.zeros=TRUE)

# energy substitution elasticities
write.gams.variable(fc,"se_en",se_en,comments="substitution elasticity - primary energy bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_ene",se_ene,comments="substitution elasticity - electricity and primary energy bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_cen",se_cen,comments="substitution elasticity - consumption primary energy bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_cene",se_cene,comments="substitution elasticity - consumption electricity and energy bundle",write.zeros=TRUE)

# labor supply elasticities
write.gams.variable(fc,"lse_sub",lse_sub,comments="labor supply elasticity - substitution",write.zeros=TRUE)
write.gams.variable(fc,"lse_inc",lse_inc,comments="labor supply elasticity - income",write.zeros=TRUE)

# household consumption elasticities
write.gams.variable(fc,"se_c",se_c,comments="substitution elasticity - consumption bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_cme",se_cme,comments="substitution elasticity - consumption materials and energy bundle",write.zeros=TRUE)
write.gams.variable(fc,"se_cm",se_cm,comments="substitution elasticity - consumption materials bundle",write.zeros=TRUE)
write.gams.variable(fc,"eta(r,h)",eta,comments="substitution elasticity - inverse of consumption over time",write.zeros=TRUE)

# fixed factor shares in resource and agriculture sectors
write.gams.variable(fc,"res_share(r,s)$sres(s)",res_share,comments="fixed factor resource share of capital for resource extraction sectors",write.zeros=TRUE)
write.gams.variable(fc,"res_share(r,s)$sagf(s)",land_share,comments="fixed factor resource share of capital for agriculture and forestry sectors",write.zeros=TRUE)

# resource extraction sector supply elasticities
write.gams.variable(fc,"se_res_sr",se_res_sr,comments="short-run supply elasticity for resoursce extraction sectors",write.zeros=TRUE)
write.gams.variable(fc,"se_res_lr",se_res_lr,comments="long-run supply elasticity for resoursce extraction sectors",write.zeros=TRUE)
write.gams.variable(fc,"rate_sr_lr",rate_sr_lr,comments="rate of convergence to long-run supply elasticity",write.zeros=TRUE)

# macroeconomic variables
write.gams.variable(fc,"gamma",gamma,comments="rate - population growth rate",write.zeros=TRUE)
write.gams.variable(fc,"omega",omega,comments="rate - labor productivity growth",write.zeros=TRUE)
write.gams.variable(fc,"delta",delta,comments="rate - depreciation rate",write.zeros=TRUE)
write.gams.variable(fc,"rbar",rbar,comments="rate - benchmark interest rate",write.zeros=TRUE)

# average energy intesity growth rate
write.gams.variable(fc,"ene_growth",ene_growth,comments="energy intesity growth rate",write.zeros=TRUE)

# average growth rate of coal fired electricity generation share
write.gams.variable(fc,"col_ele_growth",col_ele_growth,comments="growth rate of coal fired electricity generation share",write.zeros=TRUE)

# household energy intensity growth
write.gams.variable(fc,"cd_ene_growth",cd_ene_growth,comments="household energy intensity growth",write.zeros=TRUE)

# close the connection to the output file
close(fc)
